z*/7y- }D5t>7 



CASE FILE 
COPY 


AN APPLICATION OF 
A LINEAR PROGRAMING TECHNIQUE 
TO NONLINEAR MINIMAX PROBLEMS 

by James R. Schiess 

Langley Research Center 
Hampton , Va. 23665 


NATIONAL AERONAUTICS AND SPACE ADMINISTRATION • WASHINGTON, D. C • NOVEMBER 1973 



1. Report No. 

NASA TN D-7294 

2. Government Accession No. 

3. Recipient's Catalog No. 

4. Title and Subtitle 


5. Report Date 

AN APPLICATION OF A LINEAR PROGRAMING TECHNIQUE 

November 1973 

TO NONLINEAR MINIMAX PROBLEMS 

6. Performing Organization Code 

7. Author(s) 


8. Performing Organization Report No. 

James R. Schiess 


L-8789 



10. Work Unit No. 

9. Performing Organization Name and Address 


501-06-01-08 

NAbA Langley Research Center 

11. Contract or Grant No. 

Hampton, Va. 23665 





13. Type of Report and Period Covered 

12. Sponsoring Agency Name and Address 


Technical Note 

National Aeronautics and bpace Administration 

14. Sponsoring Agency Code 

Washington, D.C. 20546 



15. Supplementary Notes 


16. Abstract 


A differential correction technique for solving nonlinear minimax problems is pre- 
sented. The basis of the technique is a linear programing algorithm which solves the 
linear minimax problem. By linearizing the original nonlinear equations about a nominal 
solution, both nonlinear approximation and estimation problems using the minimax norm 
may be solved iteratively. Some consideration is also given to improving convergence 
and to the treatment of problems with more than one measured quantity. A sample prob- 
lem is treated with this technique and with the least-squares differential correction method 
to illustrate the properties of the minimax solution. The results indicate that for the sam- 
ple approximation problem, the minimax technique provides better estimates than the least- 
squares method if a sufficient amount of data is used. For the sample estimation problem, 
the minimax estimates are better if the mathematical model is incomplete. 


17. Key Words (Suggested by Author(s)) 
Linear programing 
Minimax norm 
Nonlinear estimation 
Nonlinear approximation 


18. Distribution Statement 

Unclassified ■ 

- Unlimited 


19. Security Classif. (of this report! 

20. Security Classif. (of this page) 

21. No. of Pages 

22. Price* 

Unclassified 

Unclassified 

24 

Domestic, $2.75 
Foreign, $5.25 


For sale by the National Technical Information Service, Springfield, Virginia 22151 






















AN APPLICATION OF A LINEAR PROGRAMING TECHNIQUE 
TO NONLINEAR MINIMAX PROBLEMS 

By James R. Schiess 
Langley Research Center 

SUMMARY 

A differential correction technique for solving nonlinear minimax problems is pre- ‘ 
sented. The basis of the technique is a linear programing algorithm which solves the 
linear minimax problem. By linearizing the original nonlinear equations about a nominal 
solution, both nonlinear approximation and estimation problems using the minimax norm 
may be solved iteratively. Some consideration is also given to improving convergence 
and to the treatment of problems with more than one measured quantity. A sample problem 
is treated with this technique and with the least-squares differential correction method to 
illustrate the properties of the minimax solution. The results indicate that for the sample 
approximation problem, the minimax technique provides better estimates than the least- 
squares method if a sufficient amount of data is used. For the sample estimation problem, 
the minimax estimates are better if the mathematical model is incomplete. 

INTRODUCTION 

The purpose of this paper is to present a viable approach to solving nonlinear approxi- 
mation and estimation problems by using the minimax norm. The linear approximation of, 
functions using the minimax (Chebyshev) norm has been thoroughly studied; the resulting 
approximation errors are uniformly bounded, which for many applications is highly de- 
sirable. Numerous methods for solving this problem can be found in the literature (refs. 1 
to 5). In particular, Barrodale and Young (ref. 4) treat the problem as a linear program- 
ing problem and apply a modified simplex algorithm to solve the corresponding dual 
problem (ref. 6). In this paper, the nonlinear minimax approximation problem is solved 
by linearizing the function about an assumed solution and iteratively applying the Barrodale- 
Young algorithm to improve the solution. 

Approximation problems involve the fitting of a functional form to a known continuous 
function so that a specified norm of the errors is minimized (e.g., fitting a cubic polynomial 
to sin x by using the least- squares norm). If the known function is replaced by measure- 
ments corrupted by noise and available at discrete values of the independent variable, the 
problem can be classified as an estimation problem. Therefore, by using the technique 



developed for solving nonlinear minimax approximation problems, the problem of estimating 
the parameters of a system of nonlinear equations in the presence of noisy measurements 
is also treated. 

Consideration is also given to accelerating the convergence of the technique and to 
the treatment of problems involving more than one measurement type. The accelerated 
convergence method was originally presented by Osborne and Watson (ref. 5 ) who proved 
the necessary conditions for the method to converge. Then, in order to treat problems 
» using two or more measurement types, a vector approach to minimax problems is presented. 
“ Coincident with this presentation, the manner in which data are weighted for estimation 
•problems is discussed. 

Finally, a sample problem is presented to demonstrate the application of the 
Barrodale- Young algorithm to nonlinear problems. The minimax results for both approxi- 
mation and estimation problems are presented and are compared with the results obtained 
from the least- squares differential correction method. 

SYMBOLS 


A 

AA 


F(A, t) 


n x 1 vector of unknown coefficients 
n x 1 vector of corrections to A 
jth component of A 

real-valued linear approximating function, 

L Yi (t) 

j-1 


f(t) real-valued function to be approximated; for the estimation problem, the 

uncorrupted measurement 

g. (t) jth real-valued function in linear approximating function set 

I a closed interval of the real line 

J2, J3, J4 oblateness coefficients in sample problem 
l number of measurement types at each measurement time 


2 



N 

n 

T 

t 

W 

x 0’ y 0’ z 0 
x 0’ y 0’ z 0 


a 

6 


Ht) 


number of measurements 
number of unknown coefficients 
finite subset of points in I 
independent variable 
maximum absolute error 

initial position coordinates of vehicle in sample problem, km 
initial velocity coordinates of vehicle in sample problem, km/sec 
nonnegative unknown parameters in primal problem 
noise on the ith measurement 
range rate of vehicle, km/sec 


LINEAR MINIMAX APPROXIMATION 


Before considering the nonlinear problem a brief discussion of the linear problem 
will be presented for later reference. 

Let f(t) be a known real-valued function defined on a finite subset T = jtptg, . . . 

tj^j of an interval I of the real line, and let |gj(t)j ^ be a set of n (n < N) real- 

n 

valued continuous functions defined on the same interval. Let F(A, t) = / , a-g,(t) for 

H 

r T 

the real vector A = j^a^, a2> ... , a n J . Then, the linear minimax approximation prob- 
lem consists of finding the vector A which minimizes 



Let A* be the minimizing vector and r* be the corresponding minimum value of 
expression (1). An approximation F(A*, t) that minimizes expression (1) will be called 
a best minimax approximation to f(t). 

An obvious advantageous characteristic of the minimax approximation is that once 
A* and, therefore, r* have been determined, then the error given by 

r(t.) = f(t.) - F(A*, t.) 

, at every point t. does not exceed r* in absolute value. In other words, the minimax 
norm provides a uniform approximation to the function. In contrast to this, the standard 
least- squares technique gives an approximation which fulfills a specified criterion on the 
entire set T without indicating the goodness of the approximation at any individual point. 

Under certain conditions, minimax approximations possess a characterization which 
is often used to derive methods for obtaining the approximation. First, it is necessary to 
define the Haar condition. 


Definition of Haar condition (ref. 3): A set of functions 



defined on the 


interval [a, b ] satisfy the Haar condition if each function is continuous and if, for any 
value of t e (a, b] , every set of n vectors of the following form is independent: 


[g^t), g 2 (t), . . . , g n (t)] T 


That is, the determinant 

MV MV • • • MV 
MV MV • • • MV 

MV MV • • • MV 


is nonzero for all distinct values of t^. 

The set of functions |l, x, x^, . . . , x n j, which are used in polynomial approxima- 
tions, is an example of a set fulfilling the Haar condition. The nonzero determinant for 
these functions is the well-known Vandermonde determinant. 
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The characterization of linear minimax approximations is stated in the Alternation 
Theorem as follows: 

Alternation Theorem (ref. 3): Let {g,} be a set of functions defined on fa, b] 

1 i=l - 

and satisfying the Haar condition; let T be a closed subset of [a, b]. A necessary and 
sufficient condition that 


n 

F(A, t) = a-g.(t) 

j=l 

be a best minimax approximation on T to a continuous function f(t) is that the error 
function 


r(t) = f(t) - F(A, t) 

exhibit at least n + 1 alternations on T. That is, 
r(t A ) = - r(t i+1 ) = ± r* 
for t Q < tj < . . . < t n , t. e T. 

In other words, under the hypotheses of the theorem, the error function attains the 
maximum value at least n + 1 times with alternating signs if the best approximation has 
been obtained. By using the converse of this theorem, several algorithms, such as the 
Remes single- exchange algorithm (ref. 3), construct the best minimax approximation by 
solving the n + 1 alternant equations for the coefficient set A* and the maximum error 
r*. These algorithms, however, are hampered by the fact that the proper selection of 
the points jtp, tj, . . . , t R | is unknown. Therefore, the major effort of the algorithms 
is to determine this set of points. 

It should be observed that it may not be practical to determine if the Haar condition 
holds for a particular set of functions. Because of this difficulty, other approaches for 
obtaining minimax approximations are generally more useful. 

BARRODALE -YOUNG ALGORITHM 

In the Barrodale- Young algorithm (ref. 4), the linear minimax problem is rewritten 
so that it may be treated as a linear programing problem. The algorithm is not based 
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on the Alternation Theorem or on any other special characteristics of the minimax problem. 
The dual problem to the original linear programing problem is then solved by using the 
simplex algorithm. 

First, in order to insure the nonnegativity of the parameters being sought, define the 
, , n+1 

set la. 1 by 

1 Ji j=l 


a , . = Max (o, -Min a.\ 
n+1 l lSjSn v 


°rV“n+i <1£i£n) 


Further, let 


= e(t.) = F(A, tj) - f(tj) 

and 

W = Max je^ 


By using these definitions, e . becomes 


ei = F(A, t 4 ) - f(t.) 


= H ajg] ftp - f <*i> 

i=l 


= H <*j - «n +1 > - «V 

)=1 


n n 

- H “jSj'V - ft n+l E ' f( V 
j=l 1=1 

n+1 

= Z] 

i=i 
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where 


n 

fw<v = - E ^(v 

H 

This allows the minimax problem to be written as the linear programing problem. For 
example, minimize W with respect to the nonnegative parameters where 1 < j 5 
(n + 1) and subject to the 2N constraint expressions 

n+1 

a jgj (ti)+ w > f(t 4 ) 

j =1 

n+1 

- E + w 2 - t(t i> ■ 

J 

Since in linear programing it is easier to solve a problem with many parameters 
and few constraints than one with few parameters and many constraints (ref. 6), the 
dual problem of the aforementioned is solved. 

First, the primal problem may be written in vector notation as follows: Minimize 
W = D T a 



subject to 


Ba > C 


where 

vector 


D is the (n + 2) x 1 vector D = [0, 0, . . . , 0, 1] T , a is the 
a = Oj, . . . , a n+ p W] T , B is the 2N x (n + 2) matrix, or 


(n + 2) x 1 


gjftj) . . . 

g n+l (t l> 

gl(t N ) . . . 

g n+l^W 

-g^i) . . . 

_g n+l^l^ 

~ g l^N^ “ " * 

" g n+l^W 
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and C is the 2N x 1 vector, or 

c = [%), . . . , %), -f(tj), . . . , -f(t N )] T 

In the aforementioned expressions, the symbol ">" indicates that each vector component 
fulfills this condition. 

The dual problem may immediately be written with this notation and is as follows: 
Maximize 


C T Y 


subject to 


B T Y 5 D 

where Y is a 2N x 1 vector. 

In order to solve the dual problem, nonnegative slack variables are added to the 
constraint equations to change the inequalities to equalities. A modified simplex algorithm 
is then used to solve this problem; the resulting coefficients of the slack variables in the 
solution of the dual problem are the unknown values of a. of the primal problem. 

NON UNEAR MINIMAX PROBLEMS 


The approach taken here to solve nonlinear minimax problems is the one usually 
taken in solving a nonlinear problem. That is, the nonlinear problem is linearized about 
an assumed nominal solution and an appropriate linear technique is applied. 

Let F(A, t) now be a nonlinear function of the parameter vector A. Let Aq be 
a particular "guess" of the vector which provides a nominal evaluation of F(A, t); that 
is, in some sense, F(Aq, t) approximates f(t) "closely.” Expanding F(A, t) in a 
Taylor series about Aq yields 


F(A, t) = F(A q , t) + 


— (A, t) 
dA 


A-Aq 


aA + (Higher order terms) 


where A A is the n x 1 vector, or 


aA 



Aa n] 


T 
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If the nominal F(Aq, t) is sufficiently "close" to the true solution to assume that the 
the higher order terms are small, then this series may be truncated after the linear 
term. 

By using the linearization of F(A, t), the minimax problem is to minimize 


Max If (t.) - F(A, t.)| = Max f(t.) - Jf(A q , t.) + [— (A, t.)l AA 
l<i5N 1 1 11 l£i£N 1 [ U 1 |_9A 1 J A _ A 


= Max [f(t ) - F(A 0 , t )l - \*JL (A, d AA (2) 
l<i<N L 1 U 1J id A Ja=Aq 


with respect to the difference vector A A. This is simply the linear minimax problem in 
which the function f(t^) - F(Aq, t^) is being approximated by the function 


(A, t.)l AA 

Ja=a o 


Therefore, the Barrodale-Young method can be used to obtain AA, and the nominal A Q 
can be corrected to give A^ = Aq + AA. Since the function F(A, t) is linearized, the 
vector Aj does not minimize expression (1) as desired; therefore, with the estimate of 
Aj instead of Aq, equation (2) is again minimized with respect to A A to obtain a new 
correction A A and a new estimate Ag = Aj +AA. The entire process is iterated until 
the Chebyshev norm has been sufficiently minimized. Note that, as often occurs in the 
application of the least-squares technique to a linearized problem, the process may diverge 
if the initial estimate A^ is not "close" to the best estimate of A. 

Osborne and Watson (ref. 5) have proposed a modification to this scheme and have 
proven convergence conditions for the modified scheme. In the modification, after AA 
has been obtained in the jth iteration, Osborne and Watson suggest that 

Max |f(t.) - F(A. , + /3AA, t.)| 
l<i<N 1 3-1 1 


be minimized with respect to the scalar By using the resulting value of p, the estimate 
of A becomes A^ = Aj_j +/3AA at the end of the iteration. Since p is a scalar, the 
minimization is readily obtained by using a one-dimensional technique. For the modified 
scheme, the convergence criteria are given in the following theorem: 
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Theorem (ref. 5): If the matrix VF(A) is of rank n and satisfies the Haar condi- 
tion for A in the region of interest, then the sequence Aq, Aj, Ag, . . . converges. 
Remark 1 : ^F(A) is the matrix with the N rows 


3F 

3A 


(A, tj) 


(where i = 1, 2, . . . , N) 


Remark 2 : If for a matrix of rank n every n x n submatrix is nonsingular, then the 
matrix satisfies the Haar condition. 

First, although the modified scheme is guaranteed to converge under the stated 
conditions, it should be noted that insuring that the hypotheses are true is impractical 
Specifically, for a matrix of rank n on the set of N points, proof of the Haar condition 
for one vector A requires checking the singularity of ^ submatrices. For a reason- 
able problem such as N = 40 points and n = 6 parameters, this is 3,838,380 submatrices. 
For this reason, it is easier simply to assume that the hypotheses hold. 

Second, it should be observed that the second step of the modified scheme can be 
computationally expensive, depending on the particular function F(A, t) and the tech- 
nique used to obtain p . For example, if a search technique requiring four additional 
function evaluations were applied to the sample problem presented later, this would require 
numerically integrating the differential equations four additional times during each itera- 
tion. For the current study, the unmodified algorithm without the one-dimensional search 
is considered in order to provide a valid comparison with the standard least- squares 
differential correction method. 

For nonlinear estimation problems, the preceding linearization technique can also 
be applied. The major difference between the approximation and estimation problems is 
that the function f(t) in the estimation problem is a measured quantity corrupted by 
noise. Thus, if the actual measurements are given by 


f M<V ■= 


+ e . 


where is the measurement noise and M denotes a measured quantity, then the 
problem is to minimize 


Max 
11 i!N 


[WV - F(A ' v] 


- 

1 

1 ^ 
> 

- AA 

J 

LcA 

A=An 


( 3 )- 


with respect to A A. By assuming that Aq is a good initial estimate of A and that 
F(A, t) is a perfect model of f(t), then 
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Max |f M (t.) - F(A, t.)| = Max | e.| 
l<i<N M 1 1 l<i<N 1 

In other words, if the function f(t) is perfectly modeled, the Chebyshev norm is bounded 
below by the maximum value of the measurement noise. This also indicates the major 
difficulty of minimax estimation: an extraneous data point ("wild point") drastically affects 
the norm value and consequently the estimates. Therefore, before applying a minimax 
estimation procedure, it may be advisable to edit the data judiciously. 

A limited amount of research has been devoted to the statistical aspects of minimax 
estimates. However, Howard and Kaufman (refs. 7 and 8) have derived sufficient condi- 
tions for the estimated measurements to be unbiased. In particular, if the noise samples 
e . are independent random variables with 

Max e . = - Min e . (4) 

l£i5N 1 l<ilN 1 

then the predicted measurements F(A, t) obtained from minimizing 

Max |f M (t.) - F(A, t.)| 
l<i<N M 1 1 


will be unbiased. Generally, this result is not very useful since the conditions on the 
noise given in equation (4) will not hold. Nevertheless, if the data do not contain extraneous 
points and if a large number of points are used, it is expected that the noise will approach 
the conditions in equation (4). 

FURTHER ESTIMATION CONSIDERATIONS 

Many estimation problems involve the processing of several measurement types at 
each measurement time; depending on knowledge of the measurements, it is often advisable 
or even necessary to weight the measurements in order to account for confidence in the 
data. Although the application discussions will be centered on minimax estimation based 
on unweighted scalar measurements, the more realistic problem of estimating with weighted 
vector measurements will be briefly considered in this section. 

Consider the general linear problem of estimating the n x 1 vector A based on 
a sequence of l x 1 vector measurements Yp Yg, . . . , Y m where 

Y i = GjA + e . ( i = 1, 2, . . . , m) (5) 
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In equation (5) , is an l x n matrix and e ^ is the 'Cxi measurement noise 
vector. Assuming that weight may be attached to the various components of Yp the 
weighted least-squares solution of equation (5) is obtained by minimizing with respect 
to A the function 

m m 

H £ i W £ i = L (*i - G iA) T W (Yj - GjA) (6) 

i=l i=l 

2 2 2 

In equation (6), W is usually an l x l diagonal matrix, or W = Diag (w7, w 2 , . . . , w^,), 
where w. (1 5 i < -C) is the weight associated with the ith measurement. For this 

1 rp 

reason, W can be written as the product W = V V, where V = Diag (Wp w 2 , . . . , w^,), 
and equation (6) can be written as 

m m 

£ £ Ki-v)] 1 ' [ v < Y i- G i A >] 

i=l i=l 

The original problem of equation (5) can then be restated as 

VY i = (VG i )A + Ve. (i = 1, 2, . . . , m) (7) 

In order to estimate A of equation (5) by using a "weighted minimax" approach, 
consider instead the minimax estimate of A based on equation (7). The description of a 
general linear minimax problem stated previously requires that the function being approxi- 
mated be real- valued. Since in equations (5) and (7) the function Y^ is vector-valued, a 
straightforward approach (see, for example, ref. 8) is to treat the l components of Y i 
at the m points of the data set as the l m values of a single real-valued function. If 
Y? indicates the jth component of Y^ and g| the jth row of Gp then the weighted 
minimax solution of equation (5) is obtained by minimizing 

Max | w.(Y^ - G’Ia) | 
l<i<m ] 1 1 

with respect to A, which gives the minimax solution of equation (7): 

For weighted least- squares problems, the weight w^ is usually chosen to be the 
reciprocal of the standard deviation of the noise of the corresponding measurement. Such 
a selection not only provides a relative weighting of the measurements but also means that 
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the components of equation (7) are dimensionless. For a weighted minimax approach, the 
weights may be chosen in the same fashion initially. However, since different measurement 
types are represented in equation (7), some adjustment of the weights may be required for 
a particular problem in order to insure that the various components of equation (7) are of 
the same magnitude. In other words, although the components of equation (7) are dimen- 
sionless, they may vary considerably in magnitude due to the nature of the specific measure- 
ment types; since the minimax estimates would be influenced mainly by the component of 
largest magnitude, all of the components will affect the solution only if they are all of the 
same magnitude. 

AN APPLICATION OF THE DIFFERENTIAL CORRECTION TECHNIQUE 

The problem chosen for demonstration in this study is the determination of the initial 
conditions and one of the parameters of the equations of motion of a space vehicle. The 
equations of motion, which represent the forces acting on the vehicle due to the presence 
of an oblate earth, are given by (ref. 9) 

x(t) = - + Jggjtxft)) + J 3 g 2 (x(t)} + J 4 g 3 (x(t)) (8) 

r 3 

In equation (8), x(t) and x(t) are the position and acceleration vectors, respectively, 
of the vehicle relative to a Cartesian coordinate system centered at the earth, and r is 
the magnitude of x(t). The parameters /z and J 2 , J 3 , and J 4 are, respectively, the 
gravitational constant and oblateness coefficients. The vectors g 4 (x(t)), g 2 (x(t)), and 
g 3 (x(t)) are nonlinear functions of the position. 

For this particular example, the measurement function was chosen to be the range 
rate p(t) which is a nonlinear function of the vehicle and tracking- station positions and 
velocities (ref. 9). The tracking station was assumed to be based on a rotating earth. 
Therefore, when the minimax procedure is applied to the linearized problem, equation (2) 
becomes 


Max 

l<i<N 



- Pc 




T 

] - 

— (V 

[_3P 1 J 

AP 

C 


(9) 


where c denotes values obtained from evaluations along the calculated nominal trajectory. 
The vectors P and aP are, respectively, the unknown parameters and corrections to 
the parameters. For the cases considered here, the parameters are the initial position 
(xg,yQ,ZQ), initial velocity (xQ,yQ,Zg), and the first oblateness coefficient J 2 . The 
vector of partials in expression (9) is obtained from the chain rule 
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d/o(tj) 

9x(t) 


3 p(tj) 


0 (t^) 


where the first factor on the right is a 1x6 row vector of partials of the range rate with 
respect to the position and velocity at t^ and 4 >( t) is the 6x7 matrix of partials of 
the current position and velocity with respect to the parameters. The matrix 4 > (t) is 
obtained by numerically integrating along the nominal trajectory the partials of equation (8) 
with respect to the parameters. This is accomplished by assuming that the derivatives 
are continuous in t and in the variables, so that the order of differentiation may be 
interchanged (ref. 10); for example, 


ax 

a 

dx 

_ 3 y o_ 

By 0 

.dt . 


The classical technique for solving nonlinear approximation and estimation problems 
is the least- squares differential correction method (also called modified Newton- Rap hson 
or quasilinearization). Since this is still the standard approach used on many problems, 
it was chosen to provide a comparative illustration of the results which can be expected 
from a standard technique and the minimax technique. 

Briefly described, the least- squares differential correction technique (hereinafter 
referred to simply as "least squares") uses the same linearization about an assumed 
nominal solution as that of the minimax technique; for the present example, the linearized 
problem is 


APi = P M ( V “ Pc<V (t J 

y p Jc 


AP 


Letting 


T 

• AP 

c 



( 10 ) 


(ID 


in equation (6), the least- squares criterion requires that equation (6) now be minimized 
with respect to AP. The well-known solution to this linear problem (ref. 11) is 
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AP = 


(12) 


/ 

\ 


m 


L 


BP . 


W 





The new value of P is obtained by adding AP to the initial value of P. Since the new 
P was obtained by solving the linearized problem, it does not generally minimize the 
sum of the squares of the measurement residuals. Therefore, the new value of P is 
used to obtain a second correction from equation (12). In general, the procedure must 
be iterated until a satisfactory solution is obtained. Thus, the application of the two 
linear techniques (classical least squares and linear minimax) to a nonlinear problem is 
identical except for the method (i.e., norm) used to obtain the parameter corrections. 

For the cases considered, the "true" initial conditions are given in table I. The 
"true" solution generated with these initial conditions was used for both the approximation 
and the estimation problems. The measurements (range rate) were generated every 30 
seconds along the "true" solution. For the estimation problem, normally distributed zero 
mean noise having a standard deviation of 0.1 meter per second (0.002 percent of the 
maximum range rate) was added to the measurements. For all the cases presented, the 
initial estimates of the position and velocity components were perturbed 0.1 kilometer 
and 0.1 km/sec, respectively, from the "true" values. Each technique was iterated until 
the corresponding norm of the error attained a constant minimum value (i.e., the first 
six digits were unchanged); the technique was then considered to have converged. 

By using the measurements without any noise added, the nonlinear function repre- 
sented by the range rate evaluated along the described orbit presents an excellent test 
of the outlined approach to solving a nonlinear minimax approximation problem. The first 
approximating function used in these experiments consists of the equations of motion (8) 
with J 2 , Jg, and J 4 fixed at zero; the initial conditions become the parameters available 
for adjustment. 

The results obtained by using the minimax and, as a comparison, the least- squares 
norm on the approximation problem are given in table II. The five cases presented differ 
by the number of data points used (which is equivalent to the length of the interval of the 
independent variable t). The entries in the table are the errors in the six initial conditions. 
The table suggests that, for this particular problem at least, the initial- condition errors 
are less for the minimax approximation if a sufficient number of data points (about 70 
points) are used. As an example of the function approximation errors to be encountered, 
figure 1 presents a plot of the range- rate residuals plotted against time for both norms 
applied to the 61-point case. The minimax errors all lie within the band ±2.5833 x 10" 5 
km/sec; in contrast to this, the least-squares errors generally increase, with the largest 
error in absolute value at the final point. 
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The second approximation problem consisted of fixing Jg and at zero and 
seeking the values of the initial conditions and Jg that minimize the corresponding norm. 
The initial estimate of Jg was zero. Since Jg and are three orders of magnitude 
smaller than Jg, this approximation should be very close to the original function. 

Table III presents the results of this approximation for the same five sets of data used 
previously. However, since seven parameters are required, the 41-point approximation 
diverged for both norms due to an inadequate amount of information. As in the previous 
approximation, the minimax initial condition and Jg errors are smaller than the least- 
squares errors if a sufficient amount of data is used. It should also be noted that for 
both norms, almost the same number of iterations were required on the same cases. For 
the state-only approximation problem, five iterations were required on all cases except 
in the 71-point minimax approximation which required three iterations. On the state-plus-Jg 
problem, four iterations were required by both norms except for the 51-point case for 
which both procedures used five iterations. 

Because the present problem required an extensive amount of computer time, only 
a limited number of Monte Carlo simulations were performed to test the minimax algorithm 
as an estimation procedure. The "true" trajectory and the initial estimates of the initial 
state and Jg that were used for the approximation problems were also used for the esti- 
mation problem. However, for the estimation study, ten different sets of normally distributed 
random numbers, each having a mean of zero and a standard deviation of 0.1 meter per 
second, were added to the "true" measurements to produce ten different sets of noise- 
corrupted measurements. For all of the estimation problems discussed later (e.g., 
estimation of initial state only), each of the ten sets of noisy measurements was used by 
the minimax and least-squares procedures to obtain a set of estimates. The root-mean- 
square (rms) errors of these ten sets of estimates then provided a comparison of the two 
estimation techniques. 

Since noisy measurements cause an estimation problem to be more difficult than an 
approximation problem, a larger number of data points must be used for the estimation 
problem. For the problems considered here, 61- and 81-point data sets were used. 

Neither are a sufficient number of measurements to yield accurate estimates; however, 
the estimation procedures did converge for these cases and, therefore, provide compari- 
sons between the two estimators. 

Normally, a weighted least- squares procedure is applied to an estimation problem. 

For the present problem, no weighting was used since only one measurement type was 
considered and the constant weight which would have been applied to these measurements 
can be factored out. Thus, this study reduces to a comparison of standard (unweighted) 
least-squares and minimax techniques. 
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The first estimation problem illustrates state estimation with an incomplete model. 
("Incomplete model" refers to state equations which do not completely represent all the 
dynamics present in the measurements.) The parameters Jg, Jg, and were fixed 
at zero and only the initial state was estimated. The rms errors given in the first section 
of table IV illustrate the resulting estimates. For both estimators, the estimates are 
generally improved when the number of data points is increased; however, for both data 
sets, more of the minimax estimates are better than the least-squares estimates. 

The second problem consisted of fixing the parameters at the true values and again 
estimating only the initial state. The entries in table IV show that the least- squares 
estimates are consistently better than the minimax estimates. In addition, all of the 
least-squares estimates have improved with an increase in the number of data points; 
however, most of the minimax estimates have degraded for the larger data set. 

A comparison of the results of the first two estimation problems suggests the follow- 
ing conclusions: First, for estimation problems based on an incomplete mathematical 
model, the minimax estimator may provide better estimates than the least-squares 
estimator. This result is derived from the fact that with an incomplete model most of the 
error may be caused by the model errors rather than by the measurement noise; thus, 
the effects of random noise would be negligible. Second, for problems using a complete 
model, the least- squares estimates are better than the minimax errors because the only 
error source (aside from computational errors) is the random-measurement noise; under 
these conditions the least- squares estimates are identical to the minimum variance 
estimates (ref, 11). As indicated previously, the minimax norm is bounded below by the 
maximum noise value; thus, one measurement containing noise at the three- sigma level 
has considerably more influence on the minimax estimates than on the least- squares 
estimates. This leads to the third conclusion: the minimax estimates may degrade as 
the number of measurements is increased. The degradation of the estimates is caused by 
the fact that when more measurements are processed the noise approaches the maximum 
level more closely (say, three sigma). To an extent, this degradation may be offset by 
the fact that the data interval has also increased. Thus, a trade-off occurs between the 
effects induced by increasing error bounds and longer measurement intervals. 

For the third problem, the parameters Jg and were fixed at the true values, 

Jg was initially estimated as zero, and both the state and Jg were estimated. The 
results in table IV show that the least-squares estimator again provides better estimates 
for a complete model. However, both estimators provided improved estimates when the 
data set was increased. This problem indicates the smaller data set (61 points) is approxi- 
mately the minimum number of measurements required to obtain fairly reasonable esti- 
mates of the seven parameters. For the larger data set, the estimate improvement due 
to using the 33-percent longer interval overrides any degradation of the minimax estimates 
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induced by an increased maximum noise level. Thus, problems 2 and 3 illustrate that 
minimax estimates may improve or degrade with an increase in the number of measure- 
ments, depending on the number of parameters being estimated and the number of measure- 
ments being used. 

The results of this study were obtained by using FORTRAN programs written for 
the CDC 6000 series computers. Both techniques required approximately the same number 
of memory locations; depending on the number of parameters being estimated, the minimax 
technique required 10 to 18 percent more computer time. 

CONCLUDING REMARKS 

By linearizing the original nonlinear function, nonlinear approximation and estima- 
tion problems using the minimax norm may be treated by using the Barrodale- Young 
algorithm. This approach has been applied to a simple dynamical problem. For the cases 
in which the sample problem is treated as an approximation problem, the minimax- 
parameter estimates are better than those obtained by using the least- squares technique 
when a sufficient amount of data is used. 

By adding random noise to the measurements, the sample problem has been treated 
in a limited Monte Carlo simulation as a nonlinear estimation problem. The results suggest 
that for an incomplete model the minimax estimates are better than the least- squares 
estimates; however, for a complete model the minimax estimates are worse. Because 
of this behavior, selection of one of these norms for a particular problem requires 
deciding whether accurate parameter estimates or uniformly bounded measurement errors 
is more important. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., June 29, 1973. 
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Initial conditions 
and parameters: 

Values: 

Xq, km 

7920.0 

Vn , km 

0 

Zq, km 

0 

Xg, km/sec 

0 

y Q , km/sec 

3.67204354 

Zg, km/sec 

6.11130673 

J 2 

1.0823 x 10" 3 

Jo 

-2.3 x 10" 6 

3 

J 4 

-1.8 x 10" 6 








TABLE II 

LEAST- SQUARES AND MINIMAX ERRORS FOR TWO- BODY 


APPROXIN 

LATION PI 

=tOBLEM 

Technique 

x 0’ 

km 

y 0 > 

km 

z 0’ 

km 

x 0 , 

km/ sec 

Yq> 

km/sec 

km/sec 

41 data points 

Least squares. . 
Minimax 

5.4773 

6.6781 

-3.9365 

-1.1443 

15.2882 

13.1722 

-0.011570 

-0.011453 

0.056829 

0.053413 

-0.020330 

-0.028319 

51 data points 

Least squares. . 
Minimax 

1.6789 

1.2914 

-12.0558 

-12.8184 

21.1387 

21.6812 

-0.011516 

-0.011497 

0.065399 

0.066123 

-0.035579 

-0.036330 

61 data points 

Least squares. . 
Minimax 

-0.8422 

-1.1187 

-15.5293 

-15.6605 

22.3462 

22.1912 

-0.010049 

-0.009754 

0.063877 

0.062966 

-0.035245 

-0.035125 

71 data points 

Least squares. . 
Minimax 


-14.2756 

-11.9056 

18.9720 

15.8533 

-0.007300 

-0.005729 

0.052732 

0.044871 

-0.029565 

-0.025543 

81 data points 

Least squares. . 
Minimax 

-1.7844 

-1.2839 

-8.8373 

-5.3534 

11.5633 

7.2817 

-0.003372 

-0.001341 

0.033131 

0.022690 

-0.019207 

-0.013782 
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TABLE IV 

LEAST- SQUARES AND MINIMAX ESTIMATION RESULTS 
(a) Problem 1; state estimation; incomplete model 
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.06 x 10"3 



TIME, sec 

Figure 1.- Range-rate errors for 61-point two-body approximation. 
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